library(plyr)
library(tidyverse)
library(popdemo)
library(popbio)
Load functions.
source('R/pop01_param_poun.R')
source('R/pop03_doproj.R')
source('R/pop03_doproj_stoch.R')
Each row holds a unique set of parameters that reflect diverse recovery (conservation action) and extinction (threat) scenarios.
dt <- pop01_param_poun()
The first few rows should look like this….
| species | type | first_year | a1 | a2 | a3 | a4 | b1 | b2 | b3 | b4 | c1 | c2 | c3 | c4 | d1 | d2 | d3 | d4 | akey |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Podocnemis unifilis | base | 0.0 | 0 | 0 | 0 | 10.74 | 0.0 | 0.33 | 0 | 0 | 0 | 0.17 | 0.29 | 0 | 0 | 0 | 0.11 | 0.93 | poun_base_0 |
| Podocnemis unifilis | base | 0.1 | 0 | 0 | 0 | 10.74 | 0.1 | 0.33 | 0 | 0 | 0 | 0.17 | 0.29 | 0 | 0 | 0 | 0.11 | 0.93 | poun_base_0.1 |
| Podocnemis unifilis | base | 0.2 | 0 | 0 | 0 | 10.74 | 0.2 | 0.33 | 0 | 0 | 0 | 0.17 | 0.29 | 0 | 0 | 0 | 0.11 | 0.93 | poun_base_0.2 |
nf <- 10
dt$adultF_n <- nf
# project
dout <- plyr::ddply(dt,
c("species", "type", "first_year","akey"), .fun = pop03_doproj)
dout$arun <- 1
# Model summaries
model_sum <- dout |>
group_by(species, type, first_year, lambda, gen_time) |>
summarise(fem_t0 = max(fem_t0),
fem_min = min(fem),
fem_max = max(fem)) |>
ungroup()
lambda_n <- length(unique(model_sum$lambda)) # 50
lambda_mean <- mean(model_sum$lambda) # 0.9432
lambda_sd <- sd(model_sum$lambda) # 0.1506
lambda_min <- min(model_sum$lambda) # 0.4659
lambda_max <- max(model_sum$lambda) # 1.1539
Using survival values from base model.
t0 <- 1000 # number of eggs
n_t1 <- t0 * 0.2 # 200
n_t2 <- n_t1 * 0.333333 # 66.6
n_t3 <- n_t2 * 0.333333 # 22.2
n_t4 <- n_t3 * 0.285714 # 6.3
n_t5 <- n_t4 * 0.285714 # 1.8
n_t6 <- n_t5 * 0.93 # 1.6
n_t7 <- n_t6 * 0.93 # 1.5
n_t8 <- n_t7 * 0.93 # 1.4
n_t9 <- n_t8 * 0.93 # 1.3
ceiling(t0 / n_t5) # 552 eggs to produce one adult (maturity 5 years)
## [1] 552
ceiling(t0 / n_t6) # 593 eggs to produce one adult (maturity 6 years)
## [1] 593
ceiling(t0 / n_t7) # 638 eggs to produce one adult (maturity 7 years)
## [1] 638
ceiling(t0 / n_t8) # 686 eggs to produce one adult (maturity 8 years)
## [1] 686
ceiling(t0 / n_t9) # 737 eggs to produce one adult (maturity 9 years)
## [1] 737
ceiling(552 / 20) # 28 years for a female to lay sufficient eggs (maturity 5 years)
## [1] 28
ceiling(737 / 20) # 37 years for a female to lay sufficient eggs (maturity 9 years)
## [1] 37
# Stochastic
#data frame with runs for processing
#nruns <- 100 # 100 gives same pattern as 50
nruns <- 50
dt_stoch <- dt[rep(seq_len(nrow(dt)), nruns), ]
dt_stoch$arun <- rep(1:nruns, each = nrow(dt))
# Approx 30 minutes. 1,212,000 rows. Projections quick. Summaries slow.
dout_stoch <- plyr::ddply(dt_stoch,
c("arun", "species", "type", "first_year","akey"),
.fun = pop03_doproj_stoch)
table(dout_stoch$model)
table(dout_stoch$type)
# Export for future use
saveRDS(dout_stoch, "inst/other/dout_stoch.rds")
# Combine data for plotting
dout_all <- dplyr::bind_rows(dout |> dplyr::select(arun, model, type, first_year,
akey, ayear,
lambda, gen_time,
fem, fem_t0, fem_diff, change50_flag,
change30_flag,
double_flag) |>
dplyr::mutate(lambda_lcl = NA, lambda_ucl = NA,
lambda_sd = NA, gen_sd = NA),
dout_stoch |> dplyr::select(arun, model, type, first_year,
akey, ayear,
lambda, lambda_lcl, lambda_ucl,
lambda_sd, gen_time, gen_sd,
fem, fem_t0, fem_diff, change50_flag,
change30_flag,
double_flag))
# Limit adult female number to maximum (10 x original).
dout_all[which(dout_all$fem > (nf * 10)), 'fem' ] <- (nf * 10)
# Factors in right order
dout_all$modelf <- 1
dout_all[which(dout_all$model=="Stochastic uniform") , 'modelf'] <- 2
dout_all[which(dout_all$model=="Stochastic equal") , 'modelf'] <- 3
dout_all[which(dout_all$model=="Stochastic bad x2") , 'modelf'] <- 4
dout_all[which(dout_all$model=="Stochastic bad x4") , 'modelf'] <- 5
dout_all$modelf <- as.factor(dout_all$modelf)
levels(dout_all$modelf) <- c("Deterministic", "Stochastic uniform",
"Stochastic equal", "Stochastic bad x2",
"Stochastic worst x4")
unique(dout_all$modelf)
table(dout_all$modelf)
dout_all$typef <- 1
dout_all[which(dout_all$type=="female-hunt 2.5%") , 'typef'] <- 2
dout_all[which(dout_all$type=="female-hunt 5%") , 'typef'] <- 3
dout_all[which(dout_all$type=="female-hunt 10%") , 'typef'] <- 4
dout_all[which(dout_all$type=="female-hunt 25%") , 'typef'] <- 5
dout_all[which(dout_all$type=="female-hunt 50%") , 'typef'] <- 6
dout_all$typef <- as.factor(dout_all$typef)
levels(dout_all$typef) <- c("base", "female-hunt 2.5%",
"female-hunt 5%",
"female-hunt 10%", "female-hunt 25%",
"female-hunt 50%")
table(dout_all$typef)
# first year survival
dout_all$first_yearf <- as.factor(dout_all$first_year)
fylev <- paste("first-year\nsurvival\n", seq(0, 0.9, by = 0.1), sep = "")
levels(dout_all$first_yearf) <- fylev
# Export for future use. 17/6/2024 - 1,218,060 rows 21 columns.
saveRDS(dout_all, "inst/other/dout_all.rds")
Model summaries. RedList guidelines p38 For example, upper and lower quartiles of the projected magnitude of the future reduction (i.e., reductions with 25% and 75% probability) may be considered to represent a plausible range of projected reduction Green Status Species Status: Viability A species is considered viable in a spatial unit if application of the Regional Red List Guidelines to the population in the spatial unit would result in a categorization of ‘Least Concern’ OR ‘Near Threatened and Not Declining’.
# load data
dout_all <- readRDS("inst/other/dout_all.rds")
# Make unique model ID. boot mean is same as mean (at least to 6 decimal places)
# 300 projection models.
model_ref <- dout_all |>
group_by(akey, modelf, typef, first_year, arun) |>
summarise(lambda = min(lambda),
gen_time = min(gen_time),
yc = length(unique(ayear))) |>
ungroup() |>
group_by(akey, modelf, typef, first_year) |>
summarise(count_runs = length(unique(arun)),
count_years = min(yc),
lambda_mean = mean(lambda),
lambda_min = min(lambda),
lambda_max = max(lambda),
lambda_sd = sd(lambda),
lambda_boot_lcl = Hmisc::smean.cl.boot(lambda)["Lower"],
lambda_boot_ucl = Hmisc::smean.cl.boot(lambda)["Upper"],
lambda_q25 = quantile(lambda, probs = 0.25, na.rm = TRUE),
lambda_q75 = quantile(lambda, probs = 0.75, na.rm = TRUE),
gen_mean = mean(gen_time),
gen_min = min(gen_time),
gen_max = max(gen_time),
gen_sd = sd(gen_time),
gen_boot_lcl = Hmisc::smean.cl.boot(gen_time)["Lower"],
gen_boot_ucl = Hmisc::smean.cl.boot(gen_time)["Upper"],
gen_q25 = quantile(gen_time, probs = 0.25, na.rm = TRUE),
gen_q75 = quantile(gen_time, probs = 0.75, na.rm = TRUE),
) |>
ungroup() |>
arrange(typef, first_year, modelf) |>
mutate(modelid = paste(akey, as.numeric(modelf), sep = "_"),
modelkey = row_number()) |>
relocate(modelkey, modelid)
# Estimates of generation time
gen_mean <- model_ref |>
filter(lambda_mean >= 1, first_year < 0.5) |> pull(gen_mean) |> mean()
gen_q25 <- model_ref |>
filter(lambda_mean >= 1, first_year < 0.5) |> pull(gen_mean) |>
quantile(probs = 0.25)
gen_q75 <- model_ref |>
filter(lambda_mean >= 1, first_year < 0.5) |> pull(gen_mean) |>
quantile(probs = 0.75)
gen_3 <- gen_mean * 3
gen_3_q25 <- gen_q25 * 3
gen_3_q75 <- gen_q75 * 3
# Add population changes.
model_ref2 <- model_ref |> left_join(dout_all |>
filter(ayear %in% c(ceiling(gen_3), ceiling(gen_3_q25), ceiling(gen_3_q75), 100)) |>
group_by(akey, modelf, typef, first_year, ayear, arun) |>
summarise(fem_t0 = mean(fem_t0),
femnew = mean(fem),
femnew_diff = mean(fem_diff)
) |>
ungroup() |>
group_by(akey, modelf, typef, first_year, ayear) |>
summarise(fem_t0 = mean(fem_t0),
fem = mean(femnew),
fem_sd = sd(femnew, na.rm = TRUE),
fem_diff = mean(femnew_diff),
fem_diff_sd = sd(femnew_diff, na.rm = TRUE),
fem_diff_q25 = quantile(femnew_diff, probs = 0.25, na.rm = TRUE),
fem_diff_q75 = quantile(femnew_diff, probs = 0.75, na.rm = TRUE)) |>
ungroup() |>
mutate(ayear = paste("t", ayear, sep="")) |>
pivot_wider(names_from = ayear,
values_from = c(fem, fem_sd, fem_diff,
fem_diff_sd, fem_diff_q25, fem_diff_q75)) |>
arrange(typef, first_year, modelf)
)
# Add time to change by 50% and 30%
model_ref3 <- model_ref2 |>
left_join(dout_all |>
filter(change50_flag == 1) |>
group_by(akey, modelf, typef, first_year, arun,
change50_flag) |>
summarise(acount = n(),
change50_yf = min(ayear),
change50_yl = max(ayear)) |>
ungroup() |>
group_by(akey, modelf, typef, first_year) |>
summarise(change50_ymean = floor(mean(change50_yf)),
change50_ysd = sd(change50_yf),
change50_yq25 = floor(quantile(change50_yf, probs = 0.25, na.rm = TRUE)),
change50_yq75 = floor(quantile(change50_yf, probs = 0.75, na.rm = TRUE))) |>
ungroup()
) |> left_join(dout_all |>
filter(change30_flag == 1) |>
group_by(akey, modelf, typef, first_year, arun,
change30_flag) |>
summarise(acount = n(),
change30_yf = min(ayear),
change30_yl = max(ayear)) |>
ungroup() |>
group_by(akey, modelf, typef, first_year) |>
summarise(change30_ymean = floor(mean(change30_yf)),
change30_ysd = sd(change30_yf),
change30_yq25 = floor(quantile(change30_yf, probs = 0.25, na.rm = TRUE)),
change30_yq75 = floor(quantile(change30_yf, probs = 0.75, na.rm = TRUE))) |>
ungroup()
)
# Export for future use
saveRDS(model_ref, "inst/other/model_ref.rds")
write.csv2(model_ref, "inst/other/model_ref.csv", row.names = FALSE)
# summaries of 10050 runs.
#model_all_sum <- dout_all |>
# group_by(akey, modelf, typef, first_year, arun) |>
# summarise(count_years = length(ayear),
# count_years_unique = length(unique(ayear))) |>
# ungroup()
Plots
mean_values <- model_ref |>
group_by(typef) |>
summarise(mean_l = mean(lambda_mean))
fig_model_lambda <- model_ref |>
ggplot(aes(x = first_year, y = lambda_mean, colour = modelf)) +
geom_hline(data = mean_values, aes(yintercept = mean_l),
colour = "blue", linewidth = 0.8,
linetype = "dashed") +
geom_hline(yintercept = 1,
colour = "magenta", linewidth = 0.8,
linetype = "dashed") +
geom_point() +
stat_smooth(se = FALSE) +
scale_x_continuous("first year survival and graduation",
breaks = c(0, 0.2, 0.4, 0.6, 0.8)) +
scale_color_viridis_d() +
facet_wrap(~ typef, ncol = 6) +
theme(legend.position="top") +
guides(colour=guide_legend(nrow=2,byrow=TRUE)) +
labs(y = "lambda", colour = "model")
# save
png(file = "inst/other/fig_unifilis_lambda.png", bg = "transparent",
type = c("cairo"),
width = 9, height = 4, units = "in", res=600)
fig_model_lambda
invisible(dev.off())
Check plot.
Model lambda for Podocnemis unifilis
Projections over time.
fig_proj <- dout_all |>
ggplot(aes(x = ayear, y = fem, colour = modelf)) +
geom_point(size=0.1, alpha=0.2) +
stat_smooth(se = FALSE) +
scale_colour_viridis_d("model") +
scale_y_continuous(limits = c(0, (nf * 10))) +
labs(y="Reproductive females",
x="Time (years)",
) +
theme_bw() +
theme(plot.title.position = "plot") +
facet_grid(first_yearf ~ typef)
# save
png(file = "inst/other/fig_unifilis_projections.png", bg = "transparent",
type = c("cairo"),
width = 10, height = 8, units = "in", res=600)
fig_proj
invisible(dev.off())
Check plot.
Model scenarios for Podocnemis unifilis
Number of adult females along 80 km2 of river.
# Here use density values from "tracaja_dist_5km_4z_beforeafter.R" to establish "impact"
# Adult population before - after
an_b4 <- ceiling((80 * 1.0035150)) # before = 81 in 80 km2 of river
an_b4_lci <- ceiling((80 * 0.38838535))
an_b4_uci <- ceiling((80 * 2.5928949))
an_aft <- ceiling((80 * 0.1542282)) # after = 13 in 80 km2 of river
an_aft_lci <- ceiling((80 * 0.04021797))
an_aft_uci <- ceiling((80 * 0.5914359))